# https://chat.deepseek.com/a/chat/s/b5e3ab2d-2805-4cc5-becd-574763953b80
import numpy as np
from numpy import mat, shape, zeros, multiply, nonzero, array

class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m, 1)))
        self.b = 0.0  # 确保b是标量
        self.eCache = mat(zeros((self.m, 2)))  # 修复eCache形状

def calcEk(oS, k):
    """计算第k个样本的误差"""
    # 确保所有操作返回标量值
    fXk = float(multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k, :].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek

def selectJrand(i, m):
    """随机选择不等于i的j"""
    j = i
    while j == i:
        j = int(np.random.uniform(0, m))
    return j

def selectJ(i, oS, Ei):
    """选择第二个alpha j以获得最大步长"""
    maxK = -1; maxDeltaE = 0; Ej = 0
    # 确保赋值的是标量或适当形状的数组
    oS.eCache[i] = [1, Ei]  # 更新缓存
    #【3】找到缓存中第一列（缓存有效标志位，1表示已缓存，0表示未缓存）不是0的索引，从0开始的索引
    validEcacheList = nonzero(oS.eCache[:, 0].A)[0]  # 获取有效缓存索引
    
    if len(validEcacheList) > 1:
        for k in validEcacheList:
            if k == i: continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if deltaE > maxDeltaE:  # 选择最大步长的alpha
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:  # 【4】随机选择一个不等于i的元素j，计算误差
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
        return j, Ej

def updateEk(oS, k):
    """更新误差缓存"""
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]  # 设置有效标志

def clipAlpha(aj, H, L):
    """将alpha裁剪到[L, H]区间"""
    if aj > H:
        aj = H
    if aj < L:
        aj = L
    return aj

def innerL(i, oS):
    """SMO内循环：优化alpha[i]和alpha[j]"""
    # 【1】计算第i个样本的误差
    Ei = calcEk(oS, i)
    
    # 【2】检查alpha[i]是否违反KKT条件
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or \
       ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        
        j, Ej = selectJ(i, oS, Ei)
        alphaIold = float(oS.alphas[i])  # 转换为标量
        alphaJold = float(oS.alphas[j])  # 转换为标量
        
        # 计算alpha[j]的边界
        if oS.labelMat[i] != oS.labelMat[j]:
            L = max(0, float(oS.alphas[j]) - float(oS.alphas[i]))
            H = min(oS.C, oS.C + float(oS.alphas[j]) - float(oS.alphas[i]))
        else:
            L = max(0, float(oS.alphas[j]) + float(oS.alphas[i]) - oS.C)
            H = min(oS.C, float(oS.alphas[j]) + float(oS.alphas[i]))
        
        if L == H: 
            print("L=H")
            return 0
        
        # 计算eta（核函数计算）并确保结果为标量
        # 使用float()确保结果为标量值
        eta = float(2.0 * oS.X[i, :] * oS.X[j, :].T - 
                   oS.X[i, :] * oS.X[i, :].T - 
                   oS.X[j, :] * oS.X[j, :].T)
        
        if eta >= 0: 
            print("eta>=0")
            return 0
        
        # 更新alpha[j]
        oS.alphas[j] -= float(oS.labelMat[j]) * (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(float(oS.alphas[j]), H, L)  # 确保使用标量
        updateEk(oS, j)
        
        # 检查alpha[j]变化量
        if abs(float(oS.alphas[j]) - alphaJold) < 0.00001:
            print("j not moving enough")
            return 0
        
        # 更新alpha[i]
        oS.alphas[i] += float(oS.labelMat[j]) * float(oS.labelMat[i]) * (
            alphaJold - float(oS.alphas[j]))
        updateEk(oS, i)
        
        # 计算新的偏置项b - 确保所有操作返回标量
        # 使用float()转换所有中间结果为标量
        b1 = float(oS.b) - Ei - \
              float(oS.labelMat[i]) * (float(oS.alphas[i]) - alphaIold) * \
              float(oS.X[i, :] * oS.X[i, :].T) - \
              float(oS.labelMat[j]) * (float(oS.alphas[j]) - alphaJold) * \
              float(oS.X[i, :] * oS.X[j, :].T)
        
        b2 = float(oS.b) - Ej - \
              float(oS.labelMat[i]) * (float(oS.alphas[i]) - alphaIold) * \
              float(oS.X[i, :] * oS.X[j, :].T) - \
              float(oS.labelMat[j]) * (float(oS.alphas[j]) - alphaJold) * \
              float(oS.X[j, :] * oS.X[j, :].T)
        
        # 根据alpha值选择b
        if (0 < float(oS.alphas[i])) and (oS.C > float(oS.alphas[i])):
            oS.b = b1
        elif (0 < float(oS.alphas[j])) and (oS.C > float(oS.alphas[j])):
            oS.b = b2
        else:
            oS.b = (b1 + b2) / 2.0
        
        return 1
    else:
        return 0

def smoP(dataMatIn, classLabels, C, toler, maxIter):
    """Platt SMO外循环"""
    # 【0】确保输入数据为矩阵dataMatIn是样本数*特征矩阵，classLabels是样本数*1矩阵
    dataMatIn = mat(dataMatIn)
    classLabels = mat(classLabels).transpose()
    
    # 初始化数据结构
    oS = optStruct(dataMatIn, classLabels, C, toler)
    
    iter_count = 0
    entireSet = True
    alphaPairsChanged = 0
    
    # 主循环
    while (iter_count < maxIter) and ((alphaPairsChanged > 0) or entireSet):
        alphaPairsChanged = 0
        
        if entireSet:  # 遍历所有样本
            for i in range(oS.m):
                alphaPairsChanged += innerL(i, oS)
                print(f"fullSet, iter: {iter_count} i:{i}, pairs changed {alphaPairsChanged}")
            iter_count += 1
        else:  # 只遍历非边界样本
            # 确保使用正确的数组形状
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i, oS)
                print(f"non-bound, iter: {iter_count} i:{i}, pairs changed {alphaPairsChanged}")
            iter_count += 1
        
        # 切换模式
        if entireSet:
            entireSet = False
        elif alphaPairsChanged == 0:
            entireSet = True
    
    return oS.b, oS.alphas

def calcWs(alphas, dataArr, classLabels):
    """计算权重向量w"""
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m, n = shape(X)
    w = zeros((n, 1))
    # 确保使用正确的数组形状和类型
    alphas = array(alphas).flatten()
    labelMat = array(labelMat).flatten()
    
    for i in range(m):
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w

# 创建简单数据集
def loadSimpleDataSet():
    dataMat = np.array([[1.0, 2.1], [2.0, 1.1], [1.3, 1.0], [1.0, 1.0], [2.0, 1.0]])
    labelMat = np.array([1.0, 1.0, -1.0, -1.0, 1.0])
    return dataMat, labelMat

# 训练SVM
dataArr, labelArr = loadSimpleDataSet()



b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)

# 计算权重向量
w = calcWs(alphas, dataArr, labelArr)

print("支持向量机训练结果:")
print(f"偏置项 b: {float(b):.4f}")  # 确保b转换为浮点数
print(f"权重向量 w: [{float(w[0,0]):.4f}, {float(w[1,0]):.4f}]")

# 显示支持向量
print("\n支持向量:")
alphas_array = np.array(alphas).flatten()  # 将矩阵转换为展平的数组
for i in range(len(alphas_array)):
    if alphas_array[i] > 0:
        print(f"样本 {i}: {dataArr[i]}, alpha: {alphas_array[i]:.6f}")